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SECTION  1 
INTRODUCTION 


When  a  strategic  bomber  must  search  for  a  relocatable  target  among  a  number  of  sites, 
the  mission  planning  must  include  choosing  a  search  route.  This  routing  problem  is 
closely  relat^  to  the  classical  traveling  salesman  problem  (TSP).  However,  the  bomber 
need  not  overfly  each  site,  but  must  only  come  within  sensor  range  of  it.  An  overflight  of 
one  site  may  provide  adequate  sensor  coverage  of  nearby  sites  as  well.  Thus,  the  search 
route  planning  problem  is  better  modeled  by  the  covering  salesman  problem  (CSP) 
formulated  by  Current  and  Schilling  [1].  The  CSP  is  a  generalization  of  the  traveling 
salesman  problem  in  which  the  tour  need  not  visit  every  city,  but  every  city  must 
be  within  a  certain  distance  of  some  city  on  the  tour.  A  previous  report  presented 
a  heuristic  based  on  spacefilling  curves  for  this  problem  [2].  The  result  is  a  two- 
dimensional  route  which  passes  within  a  specified  distance  of  all  the  sites  of  interest, 
while  approximately  minimizing  the  total  distance  traveled. 

When  the  number  of  potential  sites  is  large,  the  bomber’s  time  and  fuel  constraints 
may  preclude  visiting  all  of  them.  The  planner  must  then  decide  which  of  them  are 
to  be  visited.  This  report  extends  the  above  procedure  to  include  the  selection  of  an 
appropriate  set  of  sites.  The  new  procedure  also  evaluates  alternative  routing  patterns 
and  accounts  for  the  finite  turning  radius  of  the  aircraft. 

Selection  of  the  sites  to  visit  is  handled  by  imbedding  the  heuristic  router  within  an 
optimization  program.  The  outer  program  selects  an  appropriate  set  of  sites  to  be 
visited.  The  heuristic  router  finds  a  search  route  covering  those  sites.  The  outer 
program  computes  a  score  for  the  route.  It  then  adjusts  the  sites  to  be  visited  in  an 
attempt  to  optimize  the  score  (while  satisfying  time  and  fuel  constraints)  and  repeats. 
For  this  combinatorial  optimization  problem,  a  genetic  algorithm  is  used. 

Section  2  shows  how  a  flight  path  which  accounts  for  the  finite  turning  radius  of  the 
bomber  can  be  constructed.  Section  3  describes  extensions  of  the  fast  routing  algorithm 
to  include  the  “strips”  and  “arbitrary  insertion”  heuristics.  Section  4  describes  the 
genetic  algorithm  and  its  application  to  the  search  routing  problem.  Section  5  displays 
typical  results.  A  concluding  discussion  appears  in  section  6. 
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SECTION  2 
FLIGHT  PATHS 


A  flight  plan  includes  a  set  of  ujaypoints  which  a  plane  flies  through  in  succession.  In 
this  work,  waypoints  are  two-dimensional.  When  the  waypoints  are  well  separated,  the 
flight  path  can  be  approximated  by  straight  lines  between  the  waypoints.  However,  for 
closely  spaced  waypoints,  high  aircraft  speeds,  and/or  restricted  aircraft  maneuverabil¬ 
ity,  the  finite  turn  radius  of  the  aircraft  must  be  taken  into  account. 

By  convention,  a  flight  path  consists  of  straight  lines  and  circular  arcs.  One  planning 
objective  might  be  to  find  the  shortest  such  path  whic‘.  passed  through  the  specified 
sites.  It  might  be  constructed  using  the  following  manual  procedure:  Start  with  a  flat 
map  of  the  flight  area.  Insert  a  pin  at  each  of  the  points  to  be  visited.  Drop  over  each 
pin  a  hollow  cylinder  of  a  size  corresponding  to  the  desired  turn  circle.  Thread  a  string 
around  the  cylinders  in  the  chosen  order  and  draw  it  tight,  forcing  each  cylinder  against 
its  pin.  (We  assume  the  cylinders  do  not  interfere  with  each  other.)  The  string  then 
follows  the  shortest  path  for  the  chosen  visitation  order. 

This  manual  procedure  generates  a  flight  path  which  visits  each  waypoint  halfway 
through  a  turn.  This  doubles  the  number  of  significant  locations  the  navigator  must 
keep  track  of  (search  sites  plus  turn  points).  It  requires  accurate  flying,  since  an  incor¬ 
rect  airspeed  or  bank  angle  will  lead  to  an  incorrect  turn  radius  and  possibly  a  missed 
search  site.  It  also  means  that  the  aircraft  will  never  be  level  when  overflying  the  site. 
For  these  reasons,  as  well  as  convenience,  we  have  instead  implemented  the  simpler 
waypoint  turn  algorithm:  The  aircraft  initiates  each  turn  as  it  passes  a  waypoint.  The 
aircraft  then  flies  along  a  turn  circle,  with  its  center  at  a  distance  R  from  the  the  way- 
point  and  at  right  angles  from  the  original  flight  path.  Ordinarily,  the  aircraft  turns 
toward  the  next  waypoint,  as  shown  in  figure  1.  However,  if  this  would  put  the  next 
waypoint  within  the  turn  circle,  the  turn  can  be  made  in  the  opposite  direction  instead, 
thus  adding  a  loop.  This  is  shown  in  figure  2.  These  two  cases  suffice  to  handle  arbi¬ 
trarily  difficult  routing  problems,  such  as  the  closely  spaced  waypoints  shown  in  figure 
3. 

Appendix  A  presents  some  mathematical  tools  which  were  used  to  automate  this  pro¬ 
cedure. 
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Figure  1.  Aircraft  Turns  Toward  Next  Waypoint 


Figure  2.  Aircraft  Turns  Away  from  Next  Waypoint 
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Figure  3.  Waypoint  Turn  Algorithm  Applied  to  Closely  Spaced  Waypoints 
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SECTION  3 

ROUTING  PATTERNS 


A  previous  report  [2]  presented  an  algorithm  based  on  spacefilling  curves  for  finding 
heiiristic  solutions  to  the  traveling  salesman  problem  (TSP).  This  section  describes  two 
other  TSP  heuristics  which  can  also  be  adapted  to  the  search  routing  problem. 

The  spacefilling  curve  is  a  continuous  function  from  the  unit  interval  onto  the  imit 
square.  The  TSP  heuristic  required  a  function  f{x,y),  a  pseudoinverse  of  the  space¬ 
filling  curve,  which  was  used  to  order  the  cities  to  be  visited.  Different  spacefilling 
curves  lead  to  different  fimctions.  The  function  for  the  Hilbert  curve  is  shown  in  figure 
4.  This  function  preserves  nearness,  in  that  only  nearby  points  in  the  plane  are  given 
similar  values. 

The  “strips”  method  described  by  Karp  and  Steele  [3]  is  aui  older  heuristic  for  the  planar 
TSP.  For  N  cities,  it  consists  of  the  following: 


Strips  Heuristic. 

(i)  Divide  the  area  into  \/N  horizontal  strips. 

(ii)  Sort  the  cities  along  each  strip 

(iii)  Visit  the  cities  in  the  first  strip  from  right  to  left,  then  those  in  the 
second  strip  from  left  to  right,  etc. 


The  “strips”  heuristic  can  be  implemented  in  the  same  way  as  the  spacefilling  curve 
heuristics,  by  defining  an  ordering  function  similar  to  that  shown  in  figure  5  (for  three 
strips).  For  large  numbers  of  sites  distributed  uniformly  in  the  unit  square,  the  strips 
heuristic  generates  TSP  solutions  about  22  percent  longer  than  optimal  (vs.  25  percent 
over  optimal  for  the  spacefilling  curve  heuristic  using  Sieripinski’s  curve).  The  strips 
heuristic  is  more  sensitive  to  local  variations  in  site  density. 

The  “arbitrary  insertion”  heuristic  for  the  TSP  described  by  Golden  [4]  has  also  be«i 
adapted  to  the  search  routing  problem  as  follows: 

Arbitrary  Insertion  Heuristic. 

(i)  Start  with  a  path  directly  from  the  start  to  the  finish. 

(ii)  Repeat  until  all  cities  have  been  added: 
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Figure  5.  Ordering  Function  for  Strips  Heuristic 


(a)  Arbitrarily  select  city  k  not  in  the  path  to  enter  the 
path 

(b)  Find  the  edge  (t,j)  in  the  path  which  minimizes  the 
cost  of  inserting  city  k  and  insert  it  there. 


The  arbitrary  insertion  heuristic  requires  O(N^)  steps,  while  the  spacefilling  and  strips 
heuristics  require  only  0(N  In  N)  steps. 

The  arbitrary  insertion  method  could  be  used  in  a  semiautomated  router  in  which  an 
operator  specifies  part  of  the  route  (initial  and  final  points,  and  some  of  the  intermediate 
sites)  and  the  computer  zwlds  the  remaining  sites.  The  initial  route  could  include 
multiple  visits  to  some  sites. 

The  spacefilling  and  strips  heuristics  require  a  preprocessing  step:  the  city  locations 
must  be  rotated  and  scaled  into  the  unit  square.  This  is  not  necessary  for  the  arbitrary 
insertion  heuristic. 
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SECTION  4 

GENETIC  ALGORITHMS 


The  traditional  optimization  procedures  are  calculus- based  methods  (solve  for  locations 
where  the  gradient  vanishes,  or  hill  climbing:  move  in  the  direction  of  the  local  gradi¬ 
ent),  random  searches  (generate  points  at  random  and  keep  the  best),  and  enumeration 
(generate  all  possible  points  and  keep  the  best).  Calculus-based  methods  can  be  applied 
only  to  differentiable  objective  functions,  can  be  unstable,  and  are  susceptible  to  the 
foothill  problem  -  getting  trapped  at  a  local  optimum.  The  random  and  enumerative 
methods  can  be  hopelessly  inefficient  when  the  search  space  is  large. 

A  genetic  algorithm  (GA)  is  an  alternative  optimization  procedure  inspired  by  Dar¬ 
winian  evolution.  A  conventional  GA  uses  a  bit  string  which  encodes  the  problem 
parameters.  The  GA  keeps  a  population  of  '^100  such  strings  and  applies  genetic  oper¬ 
ators  to  them  to  generate  new  strings.  Each  new  string  is  evaluated,  and  better  strings 
are  preferentially  retained  for  the  next  generation. 

GAs  are  more  generally  applicable  and  robust  than  calculus-based  methods,  and  can 
be  much  more  efficient  than  random  search  and  enumerative  methods.  GAs  have 
been  successfully  applied  to  pipeline  control  [5],  structural  optimization  [6],  recursive 
adaptive  filter  design  [7],  VLSI  circuit  layout  [8],  and  many  other  problems. 

GAs  were  developed  by  John  HoUand  and  his  colleagues  at  the  University  of  Michigan 
[9].  A  short  description  of  GAs  appears  in  the  next  two  subsections.  A  more  extensive 
discussion  may  be  found  in  reference  [10]. 


4.1  DESCRIPTION 

The  genetic  algorithms  in  this  work  follow  this  reproductive  plan: 


Genetic  Algorithm. 

(i)  Create  a  random  population  of  size  N. 

(ii)  Sort  the  population  in  ascending  order  of  fitness. 

(iii)  Repeat  until  there  are  N  -i-m  members  in  all: 

(a)  Choose  randomly  a  member  Ci,  where  1  <  i  <  N 
(uniform  distribution).  Ci  is  almost  unbiased. 
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(b)  Choose  randomly  a  member  Cj,  where  i  <  j  <  N 
(uniform  distribution).  Cj  is  biased  toward  more  fit 
members. 

(c)  Choose  randomly  a  genetic  operator  gt  (see  below). 

(d)  Create  a  new  member  C  or  C  *— 

9h{Cj)  if  9k  uses  only  one  member.  Discard  C  if 
it  is  the  same  as  either  parent. 

(e)  Determine  Stness  of  C. 

(iv)  Merge  the  new  members  into  the  population,  discarding  duplicates. 

(v)  Remove  all  but  the  last  (strongest)  N  members  of  the  population. 

(vi)  Terminate  or  go  to  step  (iii). 


m  was  set  to  .207V  to  reduce  the  sorting  effort  without  materially  delaying  entry  of  new 
members  into  the  active  population.  However,  in  a  typical  application,  determining 
the  fitness  of  new  members  dominates  the  computational  effort. 

The  most  important  genetic  operator  is  crossover,  shown  in  figure  6.  This  operator 
chooses  two  positions  along  the  bit  string  and  uses  the  enclosed  fragment  from  one 
parent  and  the  end  fragments  from  the  other  parent.  Less  important  is  mutation, 
shown  in  figure  7,  which  reverses  one  or  more  bits  chosen  at  random  from  the  single 
parent. 

The  probabilities  with  which  different  genetic  operators  are  chosen  are  themselves  ad¬ 
justed  during  the  computation.  Details  are  discussed  in  appendix  B. 


4.2  EXAMPLE 

Suppose  the  problem  is  to  find  the  value  of  x  which  maximizes 

fi2{x)  =  0^®  )  sin®(5.l7rx  -1-  0.5) 

for  the  range  0  <  x  <  1,  a  problem  suggested  by  Goldberg  and  Deb  [11]. 

In  this  case,  the  bit  string  has  a  single  field  of  16  bits.  Let  these  bits  be  the  Grav  code 
for  the  integer  i.  (The  Gray  code  is  used  instead  of  the  usual  binary  code  so  that  codes 
for  neighboring  values  will  differ  by  only  one  bit.)  x  is  then  simply  1/65536. 
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Two  crossing  points  are 
selected  at  random 

1  0  0  1  1  10  1  0  0 


o 


Figure  6.  Crossover  Operator 


A  fraction  p  of  the  bits  are 
selected  at  random  to  be  flipped 


1 _ i 


DII^^B^DE 

Hi 

1 

00011101101 


Figure  7.  Mutation  Operator 
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X 

Figure  8.  Initial  Population 


Crossover  and  mutation  are  the  only  genetic  operators  used. 

An  initial  population  of  60  is  shown  in  figure  8.  By  the  10th  generation,  shown  in  figure 
9,  members  are  clustering  around  the  local  peaks  of  the  fitness  function.  By  the  30th 
generation,  shown  in  figure  10,  all  members  are  clustered  around  the  highest  peak. 


4.3  APPLICATION  TO  ROUTING 

For  the  routing  problem,  a  4'bit  field  is  used  to  select  among  16  routing  patterns 
that  governed  the  heuristic  solution  of  the  traveling  salesman  problem.  The  patterns 
comprise  the  Sierpinski  curve,  the  Hilbert  curve,  and  14  versions  of  the  strips  heuristic 
with  different  numbers  and  orientations  of  strips. 

The  remainder  of  the  bit  string  consists  of  one  bit  per  potential  site,  which  is  set  if  the 
site  is  to  be  visited  and  zero  otherwise. 

The  “raw”  fitness  function  is  the  sum  of  the  values  of  the  visited  sites.  The  constraint 
on  flight  path  length  is  handled  by  adding  a  penalty  function  as  suggested  by  Goldberg 
[10].  Members  satisfying  the  constraint  are  said  to  be  feasible.  Details  are  discussed 
in  appendix  C. 
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Generation 


Two  genetic  operators  are  used  in  addition  to  the  crossover  and  mutation  operators 
described  above. 

The  mutate-routing  operator,  shown  in  figure  11,  replaces  the  routing  pattern  with 
another  chosen  at  random.  The  primary  motivation  for  this  is  that  changing  the  routing 
pattern  and  changing  the  visited  sites  are  significantly  different  operations,  perhaps  best 
applied  at  different  rates.  Defining  separate  mutation  operators  allows  the  operator 
probabilities  to  adapt  independently  (see  appendix  B).  Early  in  the  calculation  we  would 
expect  that  changes  in  the  routing  pattern  might  well  reduce  the  search  route  length, 
so  that  the  probability  of  selecting  the  mutate-routing  operator  should  be  moderate. 
Late  in  the  calculation  we  would  expect  that  for  a  given  member,  the  routing  pattern 
and  the  set  of  sites  to  be  visited  would  be  tuned  to  one  another.  A  change  in  routing 
pattern  would  be  unlikely  to  decrease  the  search  route  length.  However,  site  changes 
might  still  improve  the  search  (by  removing  a  site  to  shorten  the  path  length  or  adding 
a  site  to  increase  the  weight).  Therefore  it  would  be  reasonable  to  mutate  sites  more 
often  than  routing  patterns. 


Ofw  field  Is  ffltod 
with  random  bits 

I - 1 


01110100100011000 

o 

01110100100011000 

Figure  11.  Mutate-Routing  Operator 
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The  other  operator  was  motivated  by  the  observation  that  the  algorithm  frequently 
chose  routes  that  visited  relatively  low  weight  sites  while  skipping  nearby  sites  with 
higher  weights.  This  appeared  to  be  due  to  the  turn  radius  linutation.  Merely  adding 
the  higher  weight  site  may  lead  to  the  addition  of  a  loop,  whose  length  makes  the  path 
infeasible.  On  the  other  hand,  deleting  the  low  weight  site  reduces  the  value  of  the 
path.  Thus  moving  the  visit  to  the  higher  weight  site  would  require  two  mutations 
(deleting  one  site  and  adding  the  other)  but  the  intermediate  states  have  low  value  and 
would  be  unlikely  to  survive.  In  GA  terminology,  the  problem  is  deceptive. 

The  deception  is  overcome  by  combining  the  two  mutations  to  form  the  move-bit  oper¬ 
ator.  This  operator,  shown  in  figure  12,  chooses  a  visited  site  at  random  and  clears  the 
corresponding  bit  in  the  bit  string.  It  then  finds  the  preceding  and  following  set  bit 
in  the  string  and  sets  some  different  bit  in  the  interval  between  them.  This  procedure 
makes  sense  only  because  the  sites  are  rearranged  in  a  preprocessing  step  so  that  neigh¬ 
boring  bits  represent  nearby  sites.  (Actually,  the  sites  are  put  into  Sierpinski  curve 
order.) 


One  bit  Is  selectsd  at  random 


Figure  12.  Move-Bit  Operator 
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Figure  13.  Score  of  Best  Member 

The  progress  of  the  calculation  may  be  shown  by  plotting  the  raw  fitness  (sum  of  values 
of  visited  sites)  of  the  best  feasible  member  in  each  generation.  Typical  plots  are  shown 
in  figure  13.  The  fitness  increases  rapidly  at  first,  then  more  slowly.  (The  algorithm 
makes  about  80  percent  of  its  progress  toward  the  final  value  in  the  first  20  percent  of 
the  time  -  a  value  which  doesn’t  change  much  with  different  run  times.) 

One  of  the  important  tuning  parameters  is  the  population  size.  Populations  of  up 
to  several  hundred  members  were  used.  Generally,  large  populations  talce  longer  to 
converge  but  reach  better  solutions,  as  shown  in  figure  13.  However,  results  are  also 
poorer  for  very  large  populations  when  the  total  number  of  evaluations  is  held  fixed. 
This  is  due  to  the  tradeoff  between  population  variability  (i.e.,  presence  of  the  building 
blocks  for  good  solutions)  and  the  number  of  generations. 

If  the  population  is  too  small,  the  members  won’t  adequately  span  the  search  space. 
If  the  population  is  too  large,  then  (for  a  fixed  number  of  evaluations)  not  enough 
generations  will  pass  to  allow  the  best  features  of  the  various  members  to  be  combined. 
This  is  shown  by  the  results  for  a  20  site  problem  shown  in  figure  14.  Four  runs  were 
made  for  each  population  size,  with  each  limited  to  4000  evaluations.  The  figure  shows 
the  mean  score  and  a  one  standard  deviation  range.  The  optimum  population  in  this 
case  is  apparently  near  150  members.  The  results  for  a  50  site  problem  shown  in  figure 
15  suggest  that  the  optimum  population  may  increase  for  larger  problems. 
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Figure  14.  Scores  as  a  Function  of  Population  for  20  Site  Problem 
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Figure  15.  Scores  as  a  Function  of  Population  for  50  Site  Problem 
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SECTION  5 
RESULTS 


Figure  16  shows  a  sample  problem  with  20  sites.  The  path  is  required  to  start  at  the 
point  labeled  1  and  end  at  the  point  labeled  2.  The  other  sites  are  shown  as  circles 
with  areas  proportional  to  their  weight  (that  is,  the  expected  likelihood  of  finding  the 
target  there). 


Figure  16.  Sample  Problem  with  20  Sites 


Figure  17  shows  the  best  path  found  through  those  sites  with  path  length  restricted  to 
3. 

Figure  18  shows  a  second  sample  problem  with  50  sites.  The  variation  among  paths 
found  by  the  genetic  algorithm  can  be  judged  from  figures  19-21,  which  show  the  three 
best  paths  with  length  less  than  5. 
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Figure  17.  Sample  Path  Through  12  of  20  Sites 

None  of  the  above  runs  used  the  arbitrary  insertion  routing  procedure,  relying  instead  on 
spacefilling  curves  and  strips  routing  patterns.  Detailed  results  and  genetic  algorithm 
parameters  are  listed  in  table  1. 
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Figure  18.  Sample  Problem  with  50  Sites 
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Figure  21.  Path  Through  19  of  50  Sites  (Length  4.09,  Weight  26.31) 


Table  1.  Parameters  for  Sample  Results 


Cities 

Visited 

Weight 

Length 

Population 

Evaluations 

Figure 

20 

12 

11.14 

2.78 

120 

4000 

17 

50 

20 

27.23 

4.91 

200 

3000 

19 

50 

19 

26.56 

4.94 

200 

3000 

20 

50 

19 

26.31 

4.99 

200 

3000 
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SECTION  6 
DISCUSSION 


A  previous  report  presented  a  fast  heuristic  based  on  spacefilling  curves  for  the  covering 
salesman  problem.  It  finds  a  two-dimensional  route  which  passes  within  a  specified 
distance  (corresponding  to  sensor  range)  of  all  the  sites  of  interest,  while  approximately 
minimizing  the  total  distance  traveled.  This  report  extends  that  procedure  to  include 
the  selection  of  an  appropriate  subset  of  sites  when  their  number  is  large,  so  that  time 
and  fuel  constraints  preclude  visiting  all  of  them.  The  finite  turn  radius  is  also  taken 
into  account. 

This  work  disregards  sensor  range,  but  only  for  simplicity.  A  previous  report  described 
procedures  for  eliminating  some  visits  because  those  cities  are  within  sensor  range  of 
cities  which  are  visited  [2].  Those  procedures  could  be  applied  with  no  significant 
changes. 

A  bigger  effort  would  be  needed  to  permit  revisits  of  sites,  and  to  properly  accumulate 
credit  for  multiple  visits.  Neither  the  spacefilling  curve  heuristic  nor  the  strips  heuristic 
easily  allows  multiple  visits.  The  arbitrary  insertion  heuristic  could  accommodate  a 
starting  route  with  multiple  visits,  and  could  be  adapted  to  generate  such  a  route  as 
well. 

This  work  still  disregards  defenses,  terrain  avoidance,  area  searches,  waypoint  offset  (so 
the  ciircraft  does  not  pass  directly  over  a  site),  and  sensor  pointing. 

Much  of  the  computational  effort  of  this  procedure  is  in  accounting  for  the  finite  turn 
radius  of  the  aircraft.  It  might  be  improved  by  first  calculating  a  lower  bound  (the 
length  of  the  piecewise  straight  line  path)  and  discarding  the  path  if  it  is  far  from 
feasible. 
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APPENDIX  A 

FLIGHT  PATH  ALGORITHM 


A.l  BACKGROUND 

Certain  properties  of  complex  arithmetic  can  be  used  to  construct  a  flight  path  consist¬ 
ing  of  straight  lines  and  circular  arcs. 

Of  course,  addition  of  complex  numbers  corresponds  directly  to  addition  of  two- 
dimensional  vectors.  Multiplying  a  complex  number  by  a  scalar  (a  real  number)  cor¬ 
responds  to  multiplying  a  vector  by  a  scalar,  and  results  in  a  vector  with  the  same 
direction  but  in  general  a  different  length. 

Multiplying  a  complex  number  by  a  second  complex  number  corresponds  to  a  rotation 
together  with  scaling.  For  example,  multiplying  by  i  corresponds  to  a  counterclockwise 
rotation  by  90°  about  the  origin.  Multiplying  by  c*^  corresponds  to  a  counterclockwise 
rotation  by  an  angle  0  about  the  origin,  where  0  is  in  radians. 

These  concepts  can  be  used  to  construct  equations  expressing  the  desired  geometric 
relationships,  which  can  then  be  solved  for  the  unknown  values.  This  is  often  simpler 
than  using  analytic  geometry  alone. 


A.2  WAYPOINT  TURN  ALGORITHM 

Let  Si . . .  Sn  be  the  points  to  be  visited.  (In  the  following,  complex  values  are  denoted 
by  lower  case  letters,  while  upper  case  letters  denote  real  values.)  We  assume  that  no 
three  points  to  be  visited  aie  along  a  straight  line,  so  that  each  point  will  correspond 
to  a  turn.  For  now,  we  will  assume  each  turn  has  the  same  radius  R.  A  sample  flight 
path  is  shown  in  figure  22.  It  consists  of  a  straight  leg  from  bh  to  the  "turn  point”  ejt, 
a  circular  arc  centered  at  Ck  which  passes  through  the  point  sjt,  and  a  second  straight 
leg  from  bk+i  to  ejk+i. 

In  the  waypoint  turn  algorithm,  the  turn  point  coincides  with  the  visited  point: 
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Figure  22.  Flight  Path  for  Single  Turn 

This  means  the  turn  center  ct  is  at  a  distance  R  from  along  a  line  perpendicular  to 
the  preceding  straight  segment.  This  relationship  corresponds  to  the  equation 


c*  =  Cfc  +  (cjt  -  bk)j^{Ti), 


(A -2) 


where  the  upper  sign  corresponds  to  a  right  turn  and  Lk  =  |eib  —  ^ibl  is  the  length  of 
the  ibth  straight  segment. 

The  turn  direction  can  be  deduced  from  the  direction  to  the  next  site.  Specifically,  we 
expect  a  right  turn  if 


The  length  Lk+i  of  the  next  straight  segment  can  be  found  using  the  Pythagorean 
theorem, 

it+l  =  k+l-C4p-Jj"-  M-'*) 

If  the  right  hand  side  of  (A-4)  is  negative,  the  next  visited  point  is  inside  the  turn 
circle  so  that  the  flight  path  is  infeasible.  In  this  case  it  is  only  necessary  to  reverse  the 
direction  of  the  turn  at  e^  and  recalculate  cj^.  This  adds  a  loop. 
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The  equation  for  cjt  in  terms  of  the  following  straight  segment  is 


Cjb  =  h+i  +  (cfc+i  -  6ib+i)  7 - (tO- 


Solving  this  for  6jfc+i,  we  have 


h+i  = 


Cfc  ±  iek+iRfLk+i 
1  ±iR/Lii+i 


(A -5) 


(A -6) 


A.3  REMARKS 

The  turn  radius  R  has  been  assumed  constant,  but  allowing  a  different  radius  at  every 
turn  would  introduce  no  significant  difficulties. 

These  formulas  can  be  used  to  generate  flight  paths  through  specified  points  consisting 
of  straight  lines  and  circular  arcs.  To  use  the  method,  the  centers  of  the  turn  circles 
must  be  found.  A  simple  method  for  determining  these  centers,  under  the  assumption 
that  the  turn  points  coincide  with  the  points  to  be  visited,  has  been  implemented.  A 
more  elaborate  method  of  positioning  these  centers  could  generate  shorter  flight  paths. 
However,  the  difference  would  be  significant  only  for  acute  turns.  If  the  preceding  and 
following  turn  points  are  distant,  the  difference  for  a  turn  of  9  is 

6L  =  R  (sin  (0/2)  -  sin0)  =  2Rsin(0/2)(l  -  cos  (0/2))  =  0(0^)  (A  -  7) 

This  amounts  to  2R  for  a  turn  of  180®,  but  .414/?  for  a  turn  of  90®  and  only  .058/?  for 
a  turn  of  45®.  Note  that  a  good  choice  of  visitation  order  (i.e.,  a  good  solution  to  the 
traveling  salesman  problem)  will  tend  to  give  small  turn  angles. 
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APPENDIX  B 

ADAPTIVE  OPERATOR  PROBABILITIES 


The  eflBciency  of  a  genetic  algorithm  depends  on  the  settings  of  several  parameters 
including  population  size  and  the  probabilities  of  applying  each  genetic  operator.  De- 
Jong  found  robust  settings  by  using  a  test  suite  of  function  optimization  problems  [12]. 
Later,  Greffenstette  used  a  genetic  algorithm  to  find  a  better  set  of  parameter  values 
[13].  However,  these  studies  used  only  the  crossover  and  mutation  operators. 

When  additional  operators  are  used,  it  is  necessary  to  set  their  selection  probabilities. 
Setting  them  by  either  of  the  above  methods  would  require  a  very  considerable  effort. 
It  is  attractive  to  set  them  adaptively  instead. 

The  following  scheme  for  adaptively  setting  operator  probabilities  was  presented  by 
Davis  [14].  Each  new  individual  that  is  better  than  the  current  best  individual  acquires 
a  “local  delta”  equal  to  the  improvement.  At  intervals  of  I  evaluations,  each  of  the 
last  W  individuals  pass  back  to  its  parent(s)  P  times  the  sum  of  its  local  delta  and 
the  deltas  inherited  from  its  progeny  (except  that  no  delta  is  passed  back  more  than 
M  generations).  Each  individual  then  passes  its  remaining  delta  to  the  operator  which 
created  it.  The  operators  accumulate  the  deltas  in  a  vector  C.  The  vector  V  of  operator 
probabilities  is  updated  by  shifting  a  fraction  5,  so  that  V  *—  SC/\C\  +  (1  —  S)V. 

The  five  parameters  W,  /,  5,  P,  and  M  effectively  replace  the  operator  probabilities 
as  governing  parameters.  Davis  found  the  values  W  =  100,  7  =  50,  5  =  0.15,  P  =  0.15, 
and  M  =  10  to  be  effective. 

There  are  two  points  that  seem  to  make  implementation  of  this  awkward.  The  first  is 
that  even  deleted  members  would  have  to  be  kept  around,  on  the  chance  that  one  of 
their  progeny  would  pass  back  some  delta.  The  other  concerns  the  handhng  of  M,  the 
limit  on  the  number  of  generations  that  delta  is  passed  ba.ck.  If  M  could  be  disregarded, 
the  deltas  could  be  propagated  in  one  pass,  starting  with  the  most  recent  member,  with 
each  step  involving  only  the  member’s  immediate  parents.  Accounting  for  M  would  seem 
to  require  a  walk  of  the  family  tree  of  each  member  having  a  local  delta,  to  a  depth  of 
M  or  the  edge  of  the  adaptation  window  (whichever  came  first).  This  procedure  would 
seem  to  grow  exponentially  with  W/population,  and  would  be  particularly  lengthy  if 
two-parent  operators  were  successful  (leading  to  large  family  trees). 

In  addition,  one  might  want  to  reward  operators  when  they  created  individuals  that 
were  better  than  any  of  their  parents,  regardless  of  whether  they  were  better  than  the 
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previous  best  individual.  This  would  give  the  adaptation  mechanism  more  information 
to  work  with. 

Because  of  the  above  concerns,  the  following  somewhat  simpler  scheme  has  been  used 
instead. 

Each  member  has  an  associated  “responsibility  vector”  R  which  indicates  the  relative 
roll  each  genetic  operator  had  in  its  creation.  If  operator  is  used  to  create  child  c 
from  a  single  parent  p,  the  new  member  gets  a  responsibility  vector 

Rc  =  PRp^{\-  P)Uk,  {B  -  1) 

where  Uk  has  a  1  in  the  kth  position  and  Os  elsewhere.  If  the  operator  used  two  parents 
m  and  /,  the  child  gets  a  responsibility  vector 

Rc  =  P(/2m  +  Rf)l2  +  (1  -  P)Uk.  {B  -  2) 

The  operators  accumulate  credit  for  improvements  in  a  vector  C.  If  the  fitness  of  the 
new  member  exceeds  that  of  its  better  parent  by  6,  C  is  incremented  by  6Rc/\Rc\.  If 
the  new  member  is  less  fit  than  at  least  one  parent,  C  is  not  changed.  Each  generation 
(m  births),  the  vector  V  of  operator  probabilities  is  updated  by  shifting  a  fraction  S: 

y  ^5C/|C|  +  (1-5)V;  (B-3) 

and  the  credit  vector  is  attenuated: 


C  * —  fiC.  {B  —  4) 

For  this  work,  P  =  0.90,  S  =  0.30m/100,  and  p  =  0.50. 

This  scheme  follows  the  two  principles  suggested  by  Davis,  namely  that  the  probability 
of  applying  an  operator  should  be  altered  in  proportion  to  the  performance  of  the 
individuals  it  creates,  and  that  local  measures  of  performance  are  insufficient  (that  is, 
operators  creating  ancestors  of  successful  individuals  should  also  be  rewarded). 

An  example  of  operator  probability  adaptation  is  shown  in  figure  23.  Each  curve  is  an 
average  over  four  runs.  The  adaptive  procedure  confirms  the  importance  of  the  crossover 
operator,  since  its  probability  steadily  increases.  On  the  other  hand,  it  suggests  the 
“move-bit”  operator  may  be  relatively  unhelpful  (or  perhaps  helpful  only  early  in  the 
optimization  process). 
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Figure  23.  Operator  Probability  Adaptation 
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APPENDIX  C 

CONSTRAINED  OPTIMIZATION 


Suppose  we  wish  to  use  a  genetic  algorithm  to  minimize  the  function 


while  satisfying  the  constraint 


iC-1) 

gix)  <  0. 

(C7-2) 

Goldberg  suggests  a  penalty  function  approach  so  that  the  objective  function  is  replaced 
by  a  compound  objective  function: 

Mx)  =  fix)  +  a^z),  iC  -  3) 

where  z  =  g(x),  q  is  a  nonnegative  penalty  coefficient  and  $  is  a  nonnegative  penalty 
function  [10]. 

Richardson  et  al.  found  that  a  penalty  function  should  not  be  too  harsh  because  the  most 
direct  path  from  a  given  solution  to  the  optimal  solution  may  require  some  infeasible 
intermediate  structures  [15].  Therefore,  the  penalty  function  should  be  chosen  to  balance 
the  preservation  of  information  against  the  pressure  for  feasibility. 

It  is  necessary  to  penalize  infeasible  members,  but  not  so  much  that  moderately  infeasi¬ 
ble  members  cannot  contribute  information.  I  also  want  to  reweird  moderately  feasible 
members  for  satisfying  the  constraint  well.  The  penalized  fitness  function  should  have 
a  maximum  at  the  feasibility  boundary  (where  p(x)  =  0).  Thus,  I  wanted  a  penalty 
function  which  is  smooth  at  0  and  has  a  slope  there  equal  to  fig'.  On  the  other  hand, 
I  didn’t  want  members  well  inside  the  feasibility  boundary  to  be  either  penalized  or  re¬ 
warded,  so  the  penalty  function  should  be  zero  for  large  negative  values.  Members  well 
outside  the  feasibility  boundary  should  be  heavily  penalized,  so  the  penalty  function 
should  be  large  for  large  positive  values. 

I  chose  the  penalty  fimction 

$(z,u;)  =  zc*/“’,  ((7-4) 

which  is  shown  with  a  solid  line  in  figure  24. 

A  simpler  alternative  which  works  almost  as  well  is 
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^{z,w)  = 


for  z  <  0 
for  2  >  0  , 


(C-S) 


which  is  shown  with  a  dotted  line  in  figure  24.  Each  function  has  a  unit  slope  at  2  =  0,  so 
the  compound  objective  function  will  have  a  local  minimum  at  the  feasibility  boundary 
if  a  =  f/g'.  a  then  plays  the  role  of  a  Lagrange  multiplier.  Each  function  also  has  a 
second  parameter  w,  a  characteristic  width,  which  is  discussed  below. 


Figure  24.  Penalty  Functions 


Siedlecki  and  Sklansky  introduced  the  concept  of  dominance  [16].  Recall  that  we  have 
assumed  a  minimization  problem,  so  that  better  solutions  correspond  to  smaller  values 
of  /.  Solution  X,-  dominates  solution  Xj  if  /*(xi)  <  fs{xj)  for  any  positive  value  of  a. 
Evidently,  this  implies  both 

f{xi)  <  f{xj)  (C  -  6) 

and 

g(xi)  <  g{xj).  {C  -  7) 
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If  a  is  too  large  the  optimal  solution  may  be  precluded,  and  if  a  is  too  small  then 
members  will  tend  to  be  infeasible.  It  is  therefore  important  that  a  be  close  to  f  /s^.  This 
ratio  is  estimated  from  the  current  population  (specifically,  from  dominant  members  on 
either  side  of  the  feasibility  boundary)  using  a  single  state  Kalman  filter  as  follows: 


Derivative  Estimate. 

(i)  Sort  the  members  on  their  constraint  values. 

(ii)  Find  a  lit  member  near  the  feasibility  border  (the  most  lit  member 
Xf  with  constraint  —  1  <  g(xf)  <  0). 

(iii)  Find  an  infeasible  member  near  the  feasibility  border  ( the  most  lit 

member  xi  with  constraint  0  <  <1^. 

(iv)  If  the  moderately  infeasible  member  was  more  fit  ( that  is,  f(xi)  < 

f{xf)),  the  derivative  f/g'  is  estimated  as  at'  =  (f(xi)  — 
/(®/))/(l/(®«)  “  filtered  derivative  estimate  is  updated 

according  to  oi  *—  {Pa'+Ra)/{P+R),  and  its  estimated  covariance 
P  is  updated  according  to  P  *—  RPI{R  +  P)  +  Q. 

(v)  Otherwise,  P  *—  P  ■¥  Q. 


Q,  R,  and  P  are  initialized  in  the  ratio  1  : 10  : 10^^,  so  the  Kalman  filter  gain  Pf{P+R) 
is  very  high  on  the  first  iteration,  drops  quickly  thereafter,  and  doesn’t  rise  again  unless 
there  are  several  generations  in  a  row  with  no  improvements. 

For  some  runs,  the  width  w  was  kept  constant  at  2.  However,  there  was  some  evidence 
that  decreasing  w  made  the  accurate  determination  of  a  less  critical.  This  was  done  by 
initially  setting  lu  to  1,  and  thereafter  to  the  range  of  the  constraint  function  for  the 
current  population: 

W  *  Zfnax  ^tnifi’  (C*  8) 
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